      SUBROUTINE E02BAF(M,NCAP7,X,Y,W,K,WORK1,WORK2,C,SS,IFAIL)
C     NAG COPYRIGHT 1975
C     MARK 5 RELEASE
C     MARK 6 REVISED  IER-84
C     MARK 8 RE-ISSUE. IER-224 (APR 1980).
C     MARK 9A REVISED. IER-356 (NOV 1981)
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     NAG LIBRARY SUBROUTINE  E02BAF
C
C     E02BAF  COMPUTES A WEIGHTED LEAST-SQUARES APPROXIMATION
C     TO AN ARBITRARY SET OF DATA POINTS BY A CUBIC SPLINE
C     WITH KNOTS PRESCRIBED BY THE USER.  CUBIC SPLINE
C     INTERPOLATION CAN ALSO BE CARRIED OUT.
C
C     COX-DE BOOR METHOD FOR EVALUATING B-SPLINES WITH
C     ADAPTATION OF GENTLEMAN*S PLANE ROTATION SCHEME FOR
C     SOLVING OVER-DETERMINED LINEAR SYSTEMS.
C
C     USES NAG LIBRARY ROUTINE  P01AAF.
C
C     STARTED - 1973.
C     COMPLETED - 1976.
C     AUTHOR - MGC AND JGH.
C
C     REDESIGNED TO USE CLASSICAL GIVENS ROTATIONS IN
C     ORDER TO AVOID THE OCCASIONAL UNDERFLOW (AND HENCE
C     OVERFLOW) PROBLEMS EXPERIENCED BY GENTLEMAN*S 3-
C     MULTIPLICATION PLANE ROTATION SCHEME
C
C     WORK1  AND  WORK2  ARE WORKSPACE AREAS.
C     WORK1(R)  CONTAINS THE VALUE OF THE  R TH  DISTINCT DATA
C     ABSCISSA AND, SUBSEQUENTLY, FOR  R = 1, 2, 3, 4,  THE
C     VALUES OF THE NON-ZERO B-SPLINES FOR EACH SUCCESSIVE
C     ABSCISSA VALUE.
C     WORK2(L, J)  CONTAINS, FOR  L = 1, 2, 3, 4,  THE VALUE OF
C     THE  J TH  ELEMENT IN THE  L TH  DIAGONAL OF THE
C     UPPER TRIANGULAR MATRIX OF BANDWIDTH  4  IN THE
C     TRIANGULAR SYSTEM DEFINING THE B-SPLINE COEFFICIENTS.
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02BAF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  SS
      INTEGER           IFAIL, M, NCAP7
C     .. Array Arguments ..
      DOUBLE PRECISION  C(NCAP7), K(NCAP7), W(M), WORK1(M),
     *                  WORK2(4,NCAP7), X(M), Y(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  ACOL, AROW, CCOL, COSINE, CROW, D, D4, D5, D6,
     *                  D7, D8, D9, DPRIME, E2, E3, E4, E5, K0, K1, K2,
     *                  K3, K4, K5, K6, N1, N2, N3, RELEMT, S, SIGMA,
     *                  SINE, WI, XI
      INTEGER           I, IERROR, IPLUSJ, IU, J, JOLD, JPLUSL, JREV, L,
     *                  L4, LPLUS1, LPLUSU, NCAP, NCAP3, NCAPM1, R
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, SQRT
C     .. Executable Statements ..
      IERROR = 4
C     CHECK THAT THE VALUES OF  M  AND  NCAP7  ARE REASONABLE
      IF (NCAP7.LT.8 .OR. M.LT.NCAP7-4) GO TO 420
      NCAP = NCAP7 - 7
      NCAPM1 = NCAP - 1
      NCAP3 = NCAP + 3
C
C     IN ORDER TO DEFINE THE FULL B-SPLINE BASIS, AUGMENT THE
C     PRESCRIBED INTERIOR KNOTS BY KNOTS OF MULTIPLICITY FOUR
C     AT EACH END OF THE DATA RANGE.
C
      DO 20 J = 1, 4
         I = NCAP3 + J
         K(J) = X(1)
         K(I) = X(M)
   20 CONTINUE
C
C     TEST THE VALIDITY OF THE DATA.
C
C     CHECK THAT THE KNOTS ARE ORDERED AND ARE INTERIOR
C     TO THE DATA INTERVAL.
C
      IERROR = 1
      IF (K(5).LE.X(1) .OR. K(NCAP3).GE.X(M)) GO TO 420
      DO 40 J = 4, NCAP3
         IF (K(J).GT.K(J+1)) GO TO 420
   40 CONTINUE
C
C     CHECK THAT THE WEIGHTS ARE STRICTLY POSITIVE.
C
      IERROR = 2
      DO 60 I = 1, M
         IF (W(I).LE.0.0D0) GO TO 420
   60 CONTINUE
C
C     CHECK THAT THE DATA ABSCISSAE ARE ORDERED, THEN FORM THE
C     ARRAY  WORK1  FROM THE ARRAY  X.  THE ARRAY  WORK1  CONTAINS
C     THE
C     SET OF DISTINCT DATA ABSCISSAE.
C
      IERROR = 3
      WORK1(1) = X(1)
      J = 2
      DO 80 I = 2, M
         IF (X(I).LT.WORK1(J-1)) GO TO 420
         IF (X(I).EQ.WORK1(J-1)) GO TO 80
         WORK1(J) = X(I)
         J = J + 1
   80 CONTINUE
      R = J - 1
C
C     CHECK THAT THERE ARE SUFFICIENT DISTINCT DATA ABSCISSAE FOR
C     THE PRESCRIBED NUMBER OF KNOTS.
C
      IERROR = 4
      IF (R.LT.NCAP3) GO TO 420
C
C     CHECK THE FIRST  S  AND THE LAST  S  SCHOENBERG-WHITNEY
C     CONDITIONS ( S = MIN(NCAP - 1, 4) ).
C
      IERROR = 5
      DO 100 J = 1, 4
         IF (J.GE.NCAP) GO TO 160
         I = NCAP3 - J + 1
         L = R - J + 1
         IF (WORK1(J).GE.K(J+4) .OR. K(I).GE.WORK1(L)) GO TO 420
  100 CONTINUE
C
C     CHECK ALL THE REMAINING SCHOENBERG-WHITNEY CONDITIONS.
C
      IF (NCAP.LE.5) GO TO 160
      R = R - 4
      I = 4
      DO 140 J = 5, NCAPM1
         K0 = K(J+4)
         K4 = K(J)
  120    I = I + 1
         IF (WORK1(I).LE.K4) GO TO 120
         IF (I.GT.R .OR. WORK1(I).GE.K0) GO TO 420
  140 CONTINUE
C
C     INITIALISE A BAND TRIANGULAR SYSTEM (I.E. A
C     MATRIX AND A RIGHT HAND SIDE) TO ZERO. THE
C     PROCESSING OF EACH DATA POINT IN TURN RESULTS
C     IN AN UPDATING OF THIS SYSTEM. THE SUBSEQUENT
C     SOLUTION OF THE RESULTING BAND TRIANGULAR SYSTEM
C     YIELDS THE COEFFICIENTS OF THE B-SPLINES.
C
  160 DO 200 I = 1, NCAP3
         DO 180 L = 1, 4
            WORK2(L,I) = 0.0D0
  180    CONTINUE
         C(I) = 0.0D0
  200 CONTINUE
      SIGMA = 0.0D0
      J = 0
      JOLD = 0
      DO 340 I = 1, M
C
C        FOR THE DATA POINT  (X(I), Y(I))  DETERMINE AN INTERVAL
C        K(J + 3) .LE. X .LT. K(J + 4)  CONTAINING  X(I).  (IN THE
C        CASE  J + 4 .EQ. NCAP  THE SECOND EQUALITY IS RELAXED TO
C        INCLUDE
C        EQUALITY).
C
         WI = W(I)
         XI = X(I)
  220    IF (XI.LT.K(J+4) .OR. J.GT.NCAPM1) GO TO 240
         J = J + 1
         GO TO 220
  240    IF (J.EQ.JOLD) GO TO 260
C
C        SET CERTAIN CONSTANTS RELATING TO THE INTERVAL
C        K(J + 3) .LE. X .LE. K(J + 4).
C
         K1 = K(J+1)
         K2 = K(J+2)
         K3 = K(J+3)
         K4 = K(J+4)
         K5 = K(J+5)
         K6 = K(J+6)
         D4 = 1.0D0/(K4-K1)
         D5 = 1.0D0/(K5-K2)
         D6 = 1.0D0/(K6-K3)
         D7 = 1.0D0/(K4-K2)
         D8 = 1.0D0/(K5-K3)
         D9 = 1.0D0/(K4-K3)
         JOLD = J
C
C        COMPUTE AND STORE IN  WORK1(L) (L = 1, 2, 3, 4)  THE VALUES
C        OF
C        THE FOUR NORMALIZED CUBIC B-SPLINES WHICH ARE NON-ZERO AT
C        X=X(I).
C
  260    E5 = K5 - XI
         E4 = K4 - XI
         E3 = XI - K3
         E2 = XI - K2
         N1 = WI*D9
         N2 = E3*N1*D8
         N1 = E4*N1*D7
         N3 = E3*N2*D6
         N2 = (E2*N1+E5*N2)*D5
         N1 = E4*N1*D4
         WORK1(4) = E3*N3
         WORK1(3) = E2*N2 + (K6-XI)*N3
         WORK1(2) = (XI-K1)*N1 + E5*N2
         WORK1(1) = E4*N1
         CROW = Y(I)*WI
C
C        ROTATE THIS ROW INTO THE BAND TRIANGULAR SYSTEM USING PLANE
C        ROTATIONS.
C
         DO 320 LPLUS1 = 1, 4
            L = LPLUS1 - 1
            RELEMT = WORK1(LPLUS1)
            IF (RELEMT.EQ.0.0D0) GO TO 320
            JPLUSL = J + L
            L4 = 4 - L
            D = WORK2(1,JPLUSL)
            IF (ABS(RELEMT).GE.D) DPRIME = ABS(RELEMT)
     *          *SQRT(1.0D0+(D/RELEMT)**2)
            IF (ABS(RELEMT).LT.D) DPRIME = D*SQRT(1.0D0+(RELEMT/D)**2)
            WORK2(1,JPLUSL) = DPRIME
            COSINE = D/DPRIME
            SINE = RELEMT/DPRIME
            IF (L4.LT.2) GO TO 300
            DO 280 IU = 2, L4
               LPLUSU = L + IU
               ACOL = WORK2(IU,JPLUSL)
               AROW = WORK1(LPLUSU)
               WORK2(IU,JPLUSL) = COSINE*ACOL + SINE*AROW
               WORK1(LPLUSU) = COSINE*AROW - SINE*ACOL
  280       CONTINUE
  300       CCOL = C(JPLUSL)
            C(JPLUSL) = COSINE*CCOL + SINE*CROW
            CROW = COSINE*CROW - SINE*CCOL
  320    CONTINUE
         SIGMA = SIGMA + CROW**2
  340 CONTINUE
      SS = SIGMA
C
C     SOLVE THE BAND TRIANGULAR SYSTEM FOR THE B-SPLINE
C     COEFFICIENTS. IF A DIAGONAL ELEMENT IS ZERO, AND HENCE
C     THE TRIANGULAR SYSTEM IS SINGULAR, THE IMPLICATION IS
C     THAT THE SCHOENBERG-WHITNEY CONDITIONS ARE ONLY JUST
C     SATISFIED. THUS IT IS APPROPRIATE TO EXIT IN THIS
C     CASE WITH THE SAME VALUE  (IFAIL=5)  OF THE ERROR
C     INDICATOR.
C
      L = -1
      DO 400 JREV = 1, NCAP3
         J = NCAP3 - JREV + 1
         D = WORK2(1,J)
         IF (D.EQ.0.0D0) GO TO 420
         IF (L.LT.3) L = L + 1
         S = C(J)
         IF (L.EQ.0) GO TO 380
         DO 360 I = 1, L
            IPLUSJ = I + J
            S = S - WORK2(I+1,J)*C(IPLUSJ)
  360    CONTINUE
  380    C(J) = S/D
  400 CONTINUE
      IERROR = 0
  420 IF (IERROR) 440, 460, 440
  440 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
  460 IFAIL = 0
      RETURN
      END
